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Abstract 

It is known that most of the craters on the surface of the Moon were created by the collision of minor bodies of the Solar 
System. Main Belt Asteroids, which can approach the terrestrial planets as a consequence of different types of resonance, are 
actually the main responsible for this phenomenon. 

Our aim is to investigate the impact distributions on the lunar surface that low-energy dynamics can provide. As a first approx- 
imation, we exploit the hyberbolic invariant manifolds associated with the central invariant manifold around the equilibrium 
point 1/2 of the Earth - Moon system within the framework of the Circular Restricted Three - Body Problem. Taking transit 
trajectories at several energy levels, we look for orbits intersecting the surface of the Moon and we attempt to define a rela- 
tionship between longitude and latitude of arrival and lunar craters density. Then, we add the gravitational effect of the Sun by 
considering the Bicircular Restricted Four - Body Problem. 

In the former case, as main outcome we observe a more relevant bombardment at the apex of the lunar surface, and a percentage 
of impact which is almost constant and whose value depends on the Earth - Moon distance <1em assumed. In the latter, it seems 
that the Earth - Moon and Earth - Moon - Sun relative distances and the initial phase of the Sun 6q play a crucial role on the 
impact distribution. The leading side focusing becomes more and more evident as cLem decreases and there seem to exist values 
of #o more favorable to produce impacts with the Moon. Moreover, the presence of the Sun make some trajectories to collide 
with the Earth. The corresponding percentage floats between 1 and 5 %. 

As further exploration, we assume an uniform density of impact on the lunar surface, looking for the regions in the Earth - Moon 
neighbourhood these colliding trajectories have to come from. It turns out that low-energy ejecta originated from high-energy 
impacts are also responsible of the phenomenon we are considering. 
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1 Introduction 



The surface of the Moon is constellated by impact craters of various sizes, mainly generated from the collision of objects coming 
from the Main Asteroid Belt. Indeed, the Inner Solar System can be reached by such minor bodies as a consequence of different 
types of resonance |T). The intense lunar bombardment took place between 3.8 and 4 Gy ago, being at the present day the 
meteroidal flux about 10 3 lower (2). 

The cratering process is interesting for several branches of science. First of all, by comparing densities of craters on different 
surfaces it is possible to derive the relative age of the corresponding terrains (see, for instance, [3,4,5]). Roughly speaking, the 
greater the density the older the surface. Also, the geological chronology of the terrestrial planets is now becoming more and 
more accurate thanks to the space missions that provide radiometric age estimates for different regions. This is especially true if 
we take as reference case the Moon, for which a great amount of data is now available. From this kind of analysis, new insight 
on the Solar System evolution can be obtained. As further aspect, the flux of impacts offers information on the Solar System 
minor bodies population. 

The main problem in all these studies resides on the fact that the crater formation is a phenomenon not fully understood yet. 
There does not exist a predictive, quantitative model of crater formation, that is, a reliable methodology that can be applied 
to all situations. The size of the crater that forms at the end of the excavation stage depends on the asteroid's size, speed 
and composition, on the collision angle, on the material and structure of the surface in which the crater forms and on the 
surface gravity of the target [6|. The problem in the determination of the crater's dimension concerns with the poorness of the 
experimental or observational data. This difficulty is usually overcome by extrapolating beyond experimental knowledge through 
scaling laws. 

This work regards the paths that impacting asteroids might have followed. In particular, we will deal with low-energy trajecto- 
ries, first derived in the Circular Restricted Three - Body Problem (CR3BP) framework applied to the Earth - Moon system and 
then analysed also accounting for the Sun gravitational attraction by means of the Bicircular Restricted Four - Body Problem 
(BR4BP). We assume the minor bodies to have already left the Main Asteroid Belt and we consider as main entrance to the 
Earth - Moon neighbourhood the stable invariant manifold associated with the central invariant manifold corresponding to the 
Li equilibrium point. 

We will look for the distribution of impacts that such orbits can create, paying attention to the fact that the Moon is locked in a 

1 : 1 spin-orbit resonance. In particular, we wonder if, for the range of energy under consideration, the Moon acts as a shield for 
the Earth or if the greatest concentration of collisions still takes place on the leading side of the surface, as other authors have 
pointed out with different approaches. See, for example, (7][8]|9). 

In a second step, from a backward integration, we attempt at discovering any other gate that can lead to a lunar impact within 
low-energy regimes. 

We recall that due to the small values of energy we consider, the impacts obtained can yield to at most 40 km in diameter 
craters. This value has been computed by applying the scaling laws of Melosh 1 6| to the Moon's surface with an impact velocity 
corresponding to the escape lunar velocity (about 2.4 km/s). 

2 The Circular Restricted Three - Body Problem 

The Circular Restricted Three - Body Problem 1101 studies the behaviour of a particle P with infinitesimal mass moving 
under the gravitational attraction of two primaries P\ and Pi, of masses mi and mj, revolving around their center of mass in 
circular orbits. 

To remove time from the equations of motion, it is convenient to introduce a synodical reference system {O y x,y,z}, which 
rotates around the z-axis with a constant angular velocity lu equal to the mean motion n of the primaries. The origin of the 
reference frame is set at the barycenter of the system and the x-axis on the line which joins the primaries, oriented in the 
direction of the largest primary. In this way we work with mi and fixed on the x-axis. 

The units are chosen in such a way that the distance between the primaries and the modulus of the angular velocity of the 
rotating frame are unitary. We define the mass ratio /1 as fi — m Tr ^ rn2 ■ For the Earth - Moon system fi — 0.012150582. 
The equations of motion for P can be written as 



x — 2y = x 



■3 — (x - fi) 3(2: + 1 - n), 



y + 2x = y 



2 = — 




(2.1) 



where r\ — [(x — n) 2 + y 2 + z 2 ] 2 and ri — [(x + 1 — /1) 2 + y 2 + z 2 ] 2 are the distances from P to P\ and Pi, respectively. 
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C> Ci C 2 <C <d C 3 <C <C 2 d<C <c 3 



Fig. 2.1 Zero-velocity curves (the intersection of the zero-velocity surfaces with the {z = 0} plane) for /i > 0. The motion is forbidden in 
the filled areas. The tick marks on the horizontal axis show the position of the primaries: Pi on the positive x--axis and P 2 on the negative 
x-axis. The case C < C\ = C$ is not displayed since the motion is allowed everywhere. 

The system 02.1b has a first integral, the Jacobi integral, which is given by 

2,2, 2(1 — jLt) 2/1 . . /.2 .2 i -2\ /0 n 

C = x +y +— H h (1 - WjW - [x +y +z . (2.2) 

n r 2 V / 

In the synodical reference system, there exist five equilibrium or libration points. Three of them, the collinear ones, are in the 
line joining the primaries and are usually denoted by L\, L 2 and L3. If xj Ji (i = 1, 2, 3) denotes the abscissa of the three collinear 
points, we assume that 

X L 2 < fl-K x Ll < fi < x L3 . 
Let Ci be the value of the Jacobi constant at the Lj equilibrium point. We have the following relation 

C\ > C 2 > C3 > C4 = C5. 

Depending on the value of the Jacobi constant, it is possible to define where the particle can move in the configuration space. 
These regions are known as Hill's regions and their boundaries are the zero-velocity surfaces. For a given mass parameter, there 
exist five different geometric configurations, corresponding to five different energy levels. If C > Ci, the infinitesimal mass 
can just move either in a neighbourhood around the largest primary or in a small neighbourhood around the smallest one. If 
O2 < C < Ci, it can move from the neighbourhood of one primary to the neighbourhood of the other one. If C3 < C < C2, 
it occurs the so-called bottleneck configuration, that is, the accessible region opens up beyond 7712. On the other hand, P can 
go toward L\ and L3 when C4 < C < C3. Finally, if C < C4 there are not forbidden regions. In Fig. 12.11 we represent the 
intersections of the zero-velocity surfaces with the {z = 0} plane. These intersections are known as zero-velocity curves. 

The collinear libration points behave as the product of two centers by a saddle. This means that around a collinear point we 
deal with bounded orbits, which are due to the central part and also with escape trajectories, which depart exponentially from 
the neighbourhood of the collinear point for t — > ±00 and are due to the saddle component. The former kind of motion belongs 
to the central invariant manifold, the latter to the hyperbolic invariant manifolds associated with the central invariant one. The 
hyperbolic manifolds consist, in particular, in one stable and one unstable. 

To be more precise, when we consider all the energy levels, the center x center part gives rise to 4-dimensional central 
manifolds around the collinear equilibrium points 1111 . Different types of periodic and quasi-periodic orbits fill the central 
invariant manifold: we refer to them as centred orbits. On the other hand, due to the hyperbolic character, the dynamics close 
to the collinear libration points is that of an unstable equilibrium. This means that each type of central orbits around L\, Li 
and 1/3 has a stable and an unstable invariant manifold. Each manifold has two branches, a positive and a negative one. They 
look like tubes of asymptotic trajectories tending to, or departing from, the corresponding orbit. These tubes have a key role in 
the study of the natural dynamics of the libration regions. When going forwards in time, the trajectories on the stable manifold 
approach exponentially the orbit, while those on the unstable manifold depart exponentially. As a matter of fact 11211131 . these 
orbits separate two types of motion. The transit solutions are those orbits belonging to the interior of the manifold and passing 
from one region to another. The non-transit ones are those staying outside the tube and bouncing back to their departure region. 

The main idea of this work is the lunar impact trajectories to be driven by the stable component associated with each type 
of central orbit. We actually adopt a methodology that does not need to compute every family of central orbit, but for sake of 
completeness we describe the most exploited ones. Planar and vertical Lyapunov orbits, halo orbits and Lissajous orbits are all 
central orbits and underlie the approach and the results of this study. 

According to Lyapunov 's centre theorem, each centre gives risqj to a family of periodic orbits, whose period tends respec- 
tively to the frequencies related to both centers, uj p and co v , when approaching the equilibrium point. These families are known 



This statement is true unless one of frequencies is an integer multiple of the other, which happens only for a countable set of values of mass ratio. 
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Fig. 2.2 We show the Poincare section at {z = 0} corresponding to the L- 2 point of the Earth - Moon system for C = 3.142003. In 
this (x, y) projection we can distinguish clearly different types of periodic and quasi-periodic orbits, all belonging to the central invariant 
manifold associated with L 2 and thus denoted as central orbits. 



as vertical Lyapunov family and planar Lyapunov family. The quasi-periodic Lissajous orbits are those associated with two- 
dimensional tori, whose two basic frequencies u)\ and uj2 tend, respectively, to the frequencies related to both centers, u v and 
lo p , when the amplitudes tend to zero. They are characterized by an harmonic motion in the {z — 0} plane and an uncoupled 
oscillation in z-direction with different periods. Furthermore, halo orbits are 3-dimensional periodic orbits that show up at the 
first bifurcation of the planar Lyapunov family. In fact, there appear two families of halo orbits which are symmetrical with 
respect to the {z = 0} plane. They are known as north and south class halo families or also first class and second class halo 
families. In Fig. 12.21 we give a representation of the above central orbits by means of the Poincare section at {z — 0} for a given 
value of C. 

Throughout the paper, we denote the central invariant manifold corresponding to a given collinear equilibrium point Lj as 
W c (Li), and the stable and the unstable invariant manifold (the hyperbolic manifolds) associated with W c (Li) as W s (Li) and 
W u (Li), respectively. 

2.1 Transit orbits belonging to W s (I/2) 

As mentioned before, our aim is to study the role that low-energy orbits might have had in the formation of lunar impact craters. 
To this end, we assume as main channel to get to the Moon the stable invariant manifold associated with the central manifold 
around the Li point, W s (Z/2). This hypothesis is based on the fact that we admit as energy levels only those belonging to the 
third regime depicted in Fig. 12.11 In particular, we consider C3 < C < C2, that is, C £ (3.184163, 3.024150) for the mass 
parameter under study. The reason for this choice is that under either the first or the second regime, there does not exist the 
possibility that a particle coming from the Outer Solar System collides with the Moon. On the other side, by discarding the more 
energetic ones we force the asteroids to approach the Moon before arriving to the Earth. 

We need an efficient way to represent the dynamics driven by W S (Z,2) for each energy level, with no distinction on central 
orbits. The main idea we develop is to determine W s (Z/2) using only the stable invariant manifolds of the planar and vertical 
periodic orbits. In particular, for a well-defined value of C, we have 

W s (L2) C W s (Planar Lyapunov) x W s (Vertical Lyapunov) . (2.3) 

Transit trajectories of the stable invariant manifold associated with any central orbit lie inside the above product. 
In what follows, we do not explain how to compute these kinds of periodic orbits 1111 . but how we derive the initial conditions 
on the associated hyperbolic manifold. Also, we give a hint on the derivation of i2. 3b . 

2.1.1 Numerical linear approximation 

There exist different methods to compute stable and unstable manifolds associated with periodic orbits, here we consider the 
linear approximation, that makes use of the eigenvectors, corresponding to the hyperbolic directions, of the monodromy matrix 
of a given periodic orbit. The reader interested in other approaches should refer to 1 14]. 




5 



Let xo(t = 0) be the initial condition of a T-periodic orbit, <f>t the flow at time t under the CR3BP vectorfield, D(j>t its 
differential and M := _D0y(xo) the monodromy matrix. If the periodic orbit is hyperbolic, then there exist Xj, XJ 1 G Spec(M) 
such that \j G ffi/{— 1, 1}. In this case, there exist a stable and an unstable manifold, which are tangent, respectively, to the Xj 
and XJ 1 eigendirections at xo(0). 

Let vs(0) and v;y(0) be, respectively, the normalized stable and unstable eigenvector corresponding to the point xo(0) on the 
periodic orbit considered, being the normalization performed in such a way that xg j/(0) has modulus equal to 1. We recall that 
in the CR3BP case, just one hyperbolic eigenvector is sufficient to determine both branches of both manifolds, since the stable 
and the unstable directions are related by a symmetry relationship. More concretely, if (v\ , «2 , «3, i>4, V5 , vq) is the eigenvector 
associated with Xj, then («i, — V2, «3, — 1>4, «5, —vq) is the eigenvector associated with XJ . 

The linear approximation for the initial conditions of the stable and the unstable manifold at xo(0) is given, respectively, by 

x s (0) =x o (0)±ev s (0), 

Xc/ (0) =x (0)±ev c/ (0), V • 

where e is some small positive parameter. The value of e fixes the size of the displacement we are performing from the periodic 
orbit to the hyperbolic manifold, if we use the above mentioned normalization for V5(0) and vjj(0). Its value must be chosen 
in such a way to guarantee that xo(0) ± evg ^(0) are still points where the linear and nonlinear manifolds are close. However, 
it cannot be too small to prevent from rounding errors. A typical value of e = 10~ 4 (in the adimensional set of units) has been 
adopted in our computations. The sign of e determines the branch of the manifold. 
If t ^ 0, we can exploit the following relations: 

x s (t) = 0t(xo(O)) ± eD&(v S (0)), 2 - 

x^(t) = 0t(xo(O))±eD^t(v^(O)). 

The stable and unstable manifolds of the periodic orbits are 2-dimensional. Once a displacement e has been selected, given a 
point xo(0) on the periodic orbit, xgjj(t), t G [0, T], provide initial conditions on the stable/unstable manifolds. In this way, 
xo(£), t G [0, T], can be thought as one of the parameters that generate the manifolds. It is usually called the parameter along 
the orbit or phase. Sometimes it will be convenient to normalise the period so that t G [0, 2ir], The other parameter is the elapsed 
time for going, following the flow with increasing/decreasing t, from the initial condition to a certain point on the manifold. This 
time interval is usually called the parameter along the flow. 

We remark that this parametrization depends on the choice of e and on the way in which the stable/unstable direction is 
normalized. A small change in e produces an effect equivalent to a small change in xo(0), in the sense that with both changes 
we get the same orbits of the manifold. Only a small shift in the parameter along the flow will be observed. This is because the 
stable/unstable directions are transversal to the flow. 



2.1.2 Geometric behaviour 

For a fixed value of C, around a given equilibrium point there are one planar and one vertical Lyapunov periodic orbit plus 
several Lissajous orbits of different amplitudes and other types of periodic and quasi-periodic motion, whose existence depends 
on the energy level considered 1 15 1. 

Let us consider, for a well-defined energy level, the intersection between the stable invariant manifold associated with different 
central orbits and the {x — 0} plane (see Fig. 12.3b . In the (y, y) and (z, z) projections, the hyperbolic manifold associated with 
the vertical Lyapunov periodic orbit gives rise to a single closed curve, the one associated with the planar Lyapunov periodic 
orbit generates, respectively, a single closed curve and a point at the origin. On the other hand, in both projections the hyperbolic 
manifold associated with a Lissajous orbit, which has dimension 3, produces an annular region, composed by infinitely many 
closed curves chained together. Clearly, this is because a periodic orbit is a S 1 object, a Lissajous orbit is a T 2 one. If we fix the 
value of one of the two phases, say (fix, characterizing a Lissajous orbit, and let the value of the other, say <j>2, to vary in [0, 2-k] 
we get one of the closed curves forming the annular region, as shown in Fig. l2.3l on the right. 

Keeping constant the value of C, distinct Lissajous orbits are found by increasing the out-of-plane amplitude and decreasing the 
in-plane one or viceversa. The (y, y) and (z, z) projections corresponding to the hyperbolic manifolds associated with different 
Lissajous orbits with similar values of the amplitudes may intersect each other, but they tend to stay one inside the other. 
Furthermore, the projections of the hyperbolic manifold associated with the Lissajous orbit with greater out-of-plane amplitude 
will be closer to the projections associated with the vertical Lyapunov periodic orbit. In the limit case, when the Lissajous orbit 
takes vertical amplitude almost as big as that of the vertical Lyapunov periodic orbit, the two projections overlap. The same 
argument holds with respect to the planar Lyapunov periodic orbit when increasing the in-plane amplitude. This is shown in 
Fie. 1241 

If we consider other sections or other types of central orbits apart from the Lissajous ones, the same qualitative behaviour 
is found. In turn, the role of outer bound is played by the hyperbolic manifold associated with the planar Lyapunov periodic 
orbit in the (y, y) projection, by the one associated with the vertical Lyapunov periodic orbit in the (z, z) plane. This result is 
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analogous to the well-known Poincare map representation of the central manifold dynamics (see Fig. 12.2b . This explains how 
the hyperbolic manifolds associated with planar and vertical Lyapunov orbits act as energy boundaries for transit orbits lying 
inside W s/u (Li). 



2.1.3 Initial conditions corresponding to transit orbits in W s (I/2) 



Following the above considerations, we carry on the computation of transit trajectories belonging to W s (Z/2) as follows. First, 
for a fixed value of the Jacobi constant, we compute the planar and the vertical Lyapunov periodic orbit around the libration 
point, as well as the initial conditions determining the proper branch (the one that escapes from the Earth - Moon neighbourhood 
backwards in time) of their corresponding stable manifold. The initial conditions on the hyperbolic manifolds of both orbits are 
integrated backwards in time, up to their first intersection with the {x = 0} section. This procedure yields two closed curves, as 
already shown in Figs. l2.3l and l2.4l 

Next, instead of setting a grid within each closed curve, we look for an uniform distribution of random initial conditions lying 
inside them. By means of a Knuth shuffle algorithm 1 16] (see Appendix) we generate a pair of random numbers for (y, y); if the 
point is inside the closed curve in the (y, y) plane, then we generate another pair of random numbers, (z, i), that now should 
stay in the interior of the closed curve of the (z, z) projection. When both pairs fulfill the requirements, we complete the initial 
conditions setting x — and determining the x coordinate by means of the Jacobi first integral. 

We notice that taking initial conditions on the {x = 0} section means assuming the asteroids to have already left the Main 
Asteroid Belt and to have started moving in the Earth - Moon neighbourhood. Other sections apart from the {x — 0} one could 
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Fig. 2.3 (y,y) and (z,z) projections of the intersection with the {x = 0} plane of the stable manifolds associated with different types of 
orbits of energy G = 3.163, around the L2 point. On the left, on the top one can see that the projection corresponding to the manifold 
associated with the planar Lyapunov orbit is a single closed curve which contains the projections of all the other central orbits; on the 
bottom, the projection corresponding to the manifold associated with the vertical Lyapunov orbit is again a single closed curve which 
contains all the other projections. On the right, we represent the same behaviour, underlining that each closed curve constituting the 
projection of the hyperbolic manifold associated with the Lissajous orbit corresponds to a fixed value of one of the two phases. 
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have been selected. In the future, we will use the same methodology with a different choice, depending on the aspect we want 
to emphasize. For instance, the main source of colliding orbits to belong to the ecliptic plane. 



3 Impacts coming from W s (L2) 



The first exploration performed consists in numerically integrating the equations of motion of the CR3BP starting from the initial 
conditions derived as just explained. We take 10 5 points inside the (y, y) curve and, for each of these points, 10 pairs of (2, 2) 
coordinates. Hence, for each energy level, we explore the behaviour of 10 6 trajectories. This selection results in a distribution 
of initial conditions corresponding to transit orbits, which is uniform both in (y, y) and (2, 2) coordinates. However, a different 
matching among coordinates would be possible. 
Once again, we stress that the trajectories corresponding to such initial conditions will be driven by the stable component of 
W c (Z/2), without lying on W s (I/2)- They stay inside the dynamical tube generated by W s (Z/2) 1131 . 

The procedure is implemented for 10 equally spaced energy levels C in the range C3 < C < C2. In Fig. 13. II we show the Hill's 
region this energy range corresponds to, the boundary of W s (Z/2) in red and some transit trajectories in gray. 

As the more intense lunar bombardment took place some billions years ago, we must consider different values for the Earth - 
Moon distance d^M- Indeed, the Moon is receding from the Earth: the rate of recession has not been constant in the past and it 
did not behave linearly either (see, for example, I17II18II19II201 ). We take 4 values for d EM , 232400, 270400, 308400, 384400 
km, respectively. According to 1201 . they correspond approximately to 4., 3.4, 2.5 and Gy ago. 
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Fig. 2.4 (y, y) and (z, i) projections of the intersection with the {x = 0} plane of the stable manifolds associated with different types of 
orbits of energy C = 3.163, around the L 2 point. On the top, the Lissajous orbit has a considerable out-of-plane amplitude (and a very 
small in-plane one) and thus the projections of the associated stable manifold intersects the ones associated with the vertical Lyapunov 
periodic orbit. On the bottom, the Lissajous orbit has a considerable in-plane amplitude (and a very small out-of-plane one) and thus the 
projections of the associated stable manifold intersects the ones associated with the planar Lyapunov periodic orbit. 
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Fig. 3.2 We display the two ways in which a particle can escape from the Earth - Moon neighbourhood. On the left, the asteroid jumps on 
the L2 gate; on the right, it performs some loops around the Earth and then joins W U (L2) to escape. 



The maximum allowed time for impacting onto the surface of the Moon is 60 years, provided the assumption of a no longer 
life in the region under consideration. Within this time span, we get a numerical evidence that the minor bodies can behave in 
one of the following ways: 

(1) they collide with the Moon without overcoming the L\ border; 

(2) they collide with the Moon after overcoming the L\ border and thus performing several loops around the Earth; 

(3) they keep wandering around the Earth inside the area delimited by the zero-velocity surface; 

(4) they escape from the Earth - Moon neighbourhood just after jumping on the Li gate; 

(5) they exit from the Earth - Moon neighbourhood after wandering for a certain interval of time around the Earth. 

Note that just the first two cases cause the formation of craters of impact on the surface of the Moon. 
If a trajectory collides with the Moon we calculate the longitude and latitude corresponding to the site of impact, the velocity, 
the angle and the orbital elements of the osculating ellipses with respect to the Earth at the initial condition. We recall that we 
neglect the Moon's orbital inclination with respect to the Earth's orbit. 

Regarding the fourth and the fifth item (see Fig. 13.2b . the mechanism of escaping is produced by W U (L2). In the future, 
it will be interesting to analyse how the impact phenomena are fostered by homoclinic connections associated with L^. By 
finding succeeding intersections, in a given energy level, between the stable and the unstable manifold associated with the planar 
Lyapunov periodic orbit and simultaneously between the stable and the unstable manifold associated with the vertical Lyapunov 
one, it is possible to construct cycling paths, which bring the particle in and out the region demarcated by the zero-velocity 
curve. However, preliminary simulations suggest the percentage of impacts offered by this loop to be small in comparison with 
the total amount of collisions. 
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3.1 Results 

We can highlight the following outcomes: 

(a) the percentage of impacting orbits over all the initial conditions launched is 13%; 

(b) the smaller cIem, the greater the above percentage; 

(c) the amount of particles that still wander around the Earth inside the zone bounded by the zero-velocity surface after 60 years 
is 0.1%; 

(d) in all the cases of collision, the heaviest probability of impact takes place at the apex of the lunar surface (90° W, 0°). 

Point (a) and (c), in particular, reveal that the 60 years assumed are not restrictive with respect to a lunar impact. We also notice 
that in the time interval considered the most of the asteroids escapes from the region we are interested in. As already mentioned, 
it looks like just few of them are able to go back to the Earth - Moon neighbourhood later. It is reasonable to think that they 
remain in the Inner Solar System and occasionally are pushed towards the Earth again. 
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Fig. 3.3 We display the density of impact (number of impacts per unit of area normalized with respect to the total number of impacts 
obtained) computed by exploiting the CR3BP equations of motion. The surface of the Moon is discretized in squares of 15° x 15° and 4 
different values for the Earth - Moon distance are considered. On the top, d E M = 232400 km and d EM = 270400 km; on the bottom, 
<1em = 308400 km and d EM = 384400 km. The color bar indicates that the lighter the color of the square the greater the impact density. 



In Fig. 13. 31 we represent the surface of the Moon in terms of longitude and latitude. In particular, we discretize the lunar sphere 
in squares of 15° x 15°. We notice different colors associated with each square: the lighter the color the greater the number of 
collisions per unit of area normalized with respect to the total number of impacts obtained. We remark that to consider another 
discretization of the lunar surface would not bring any relevant difference from a qualitatively point of view (see Figs. 13.41 and 

1551. 

As said, we compute the orbital elements of the osculating ellipses at t = with respect to the Earth, corresponding to initial 
conditions leading to impact. This means that every set (x, y, z, x, y, z) providing a collision with the Moon is transformed to 
inertial coordinates centered at the Earth and these are then turned into orbital elements. In particular, we get the semi-major axis 
a, the eccentricity e, the inclination i, the longitude of the ascending node fl, the argument of perigee uj and the true anomaly v. 
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Fig. 3.4 We display the density of impact (number of impacts per unit of area normalized with respect to the total number of impacts 
obtained) computed by exploiting the CR3BP equations of motion. The surface of the Moon is discretized in squares of 5° x 5°. On the 
left, d E M = 232400 km; on the right, d EM = 308400 km. The color bar indicates that the lighter the color of the square the greater the 
impact density. 
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Fig. 3.5 We display the density of impact (number of impacts per unit of area normalized with respect to the total number of impacts 
obtained) computed by exploiting the CR3BP equations of motion. The surface of the Moon is discretized in squares of 10° x 10°. On the 
left, d E M = 232400 km; on the right, d EM = 308400 km. The color bar indicates that the lighter the color of the square the greater the 
impact density. 



We notice that the initial conditions are taken far enough (at least about 500000 km if we assume d^M = 384400 km) from the 
Moon to be allowed to assume a Two-Body approximation and perform this analysis. 

It turns out that the impact is more likely if (a, e, i) belong to the intervals showed in Tab. 13.11 In Fig. 13.61 we display such 
probabilities for the case d^M = 384400 km. As before, the lighter the color associated with a given (i,a)/(i,e) square, the 
greater the probability that such orbital elements would correspond to a colliding trajectory. The probability is normalized with 
respect to the total number of impacts obtained. 



Table 3.1 For each initial condition belonging to W 3 (i2) and colliding with the lunar surface, we compute the orbital elements which 
correspond to the osculating ellipse at t = with respect to the Earth. The impact is more likely if the semi-major axis a, the eccentricity e 
and the inclination i lie in the range shown here. 



a (d EM ) e 



[1.5 : 3] [0.4 : 0.7] [1.5° : 3.5°] 
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Fig. 3.6 For each initial condition belonging to W S (L2) and colliding with the lunar surface, we compute the orbital elements which 
correspond to the osculating ellipse at t = with respect to the Earth. The plot on the left shows the probability of impact considering as 
variables the inclination i and the semi-major axis a; the one on the right takes as variables i and the eccentricity e. The color bar refers to 
the number of impacts normalized with respect to the total number of impacts found. The lighter the color the greater the probability. For 
these plots the Earth - Moon distance is assumed to be d E M = 384400 km. The i and a ranges are discretized at steps of 0.5 degrees 
and d EM , respectively. The e range is discretized at steps of 0.1. 



4 Uniform density of lunar impacts: possible paths 

In the previous section, we have seen that W s (Z/2) provides a non-uniform density of impact on the surface of the Moon. The 
question that naturally arises is where minor bodies producing an uniform distribution of low-energy collisions would come 
from. 

Having this purpose in mind, for every energy level considered in Sec. [3] we create a set of initial conditions uniformly 
distributed on the lunar surface, discretized as before in terms of longitude and latitude. In this case, not only the position 
coordinates have to be well spread out, but also the velocity ones. We apply the CR3BP equations of motions to such initial 
conditions backwards in time, up to a maximum of 5 years, detecting how many trajectories arrive from the {x = 0} section 
already mentioned. 

To be more precise, if 7 € [— it, it] is a random value of latitude, ip 6 [0, 2n] a random value of longitude and pM = 
1737.53/dgM the adimensional radius of the Moon, then att = (x,y, z) are computed as: 

x = p M cos (7) cos (ip) + At - 1, y = /o M cos (7) sin (V>), z = p M sia(j). 

Both 7 and ip are generated with the Knuth shuffle algorithm |T6| (see Appendix). 

Concerning (±, y, z) at t = 0, we implement three different procedures for each 15° x 15° square in order to ensure to fulfill the 
constraint of uniform distribution on a semisphere of velocities. If g = (g x , gy,gz) = (x — p + 1, y, z) is the normal at (x, y, z) 
to the surface vector and w and h are random values in [0, 1], the three approaches can be sketched as follows. 

(1) Let f3 G [—it, it] and A £ [0, it] be random values. Then 

x = g x cos (A) cos (/3), y = g y cos (A) sin (/?), z = —g z sin (A). 

(2) Let f3 be a random value belonging to the interval [—it, it] and A = cos -1 (1 — 2w) G [0, it]. Then 

x = g x cos (A) cos (/3), y = g y cos (A) sin (/?), z = -g z sin (A). 

(3) Let 7 and ip as above, £ = 2w — 1 and 77 = 2h — 1 such that — 1 < £,r) < 1 and £ 2 + rj 2 < 1. Then 

x = (2f \/l — £ 2 — rj 2 ) cos (7) cos (ip), y — (2r]\/ 1 — £ 2 — rj 2 ) cos (7) sin (ip), z = — [1 — 2(£ 2 + rj 2 )] sin (ip). 

In any case, (x, y, z) are normalized in order to obtain the modulus of the velocity as the one satisfying the chosen C. 
In this way, we simulate the behaviour of 758640 particles for each value of C, which means 7586400 particles in total. This 
value has been chosen to be in agreement with the computations described in Sec. [3] but also to have an impact density of 
2 x 10~ 2 km~ 2 for each 15° x 15° square considered on the surface of the Moon. 
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Fig. 4.1 If the impact distribution on the surface of the Moon was uniform, initial conditions associated with this pattern inside the (y,y) 
curve wo uld c ollid e wit h the Moon. We remark that the uniform distribution would not be due only to W S (L 2 ), but also to other phenomena. 
See Fias.lT2landl4~3l 
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Fig. 4.2 We show the density of impact caused by the dynamics associated with W S (L 2 ) if the lunar craters distribution was uniform. We 
recall that the surface of the Moon is discretized in squares of 15° x 15° and that the color bar refers to the number of impacts per unit of 
area normalized with respect to the total number of impacts found. The lighter the color the greater the impact density. The Earth - Moon 
distance is assumed to be d EM = 384400 km. 



4.1 Results 



The backward simulation reveals that there exist two main dynamical channels leading to lunar collision, for the range of energy 
under study. In particular, uniform distributed impacts would come either from W s (La) or from double collision orbits with the 
surface of the Moon. 

In the first case, we notice that all the orbits getting to the {x = 0} section give rise to points which lie inside the (y, y) and 
(z, z) curves introduced in Sec. 12.1.21 This fact can be viewed as a further confirmation of the well-posed procedure adopted 
previously. With this, we mean that the strategy defined to determine W s (L2) actually represents the dynamics we were looking 
for and does not leave out any transit trajectory. The interesting point is that inside these curves we are able to note special 
patterns, that should be investigated with more detail (see Fig. 14. lb . 

The density of impact on the Moon's surface produced by W S (L2) in this case is depicted in Fig. 14.21 The reader should be 
aware that we do not expect an apex concentration as before, due to the different collocation of points inside the (y, y) and (z, z) 
projections. 

On the other hand, there exist orbits that depart from the Moon with about the lunar escape velocity and return there with the 
same speed (see Fig. 14.3b . They can travel along different paths, turning around either the Earth or the Moon one or several 
times. As explanation, we can hypothesize ejecta deriving from high-energy collisions. Such effect has already been predicted 
by other authors [21 1. 
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5 The Bicircular Restricted Four - Body Problem 

The Bicircular Restricted Four - Body Problem [22 1 considers the infinitesimal mass P to be affected by the gravitational 
attractions of three primaries. We introduce this model in order to include the effect of the Sun in the Earth - Moon system. The 
Earth and the Moon revolve in circular orbits around their common center of mass and, at the same time, this barycenter and the 
Sun move on circular orbits around the center of mass of the Earth - Moon - Sun system. 

The usual framework to deal with is the synodical reference system centered at the Earth - Moon barycenter: in this way the 
Earth and the Moon are fixed on the x-axis as before and the Sun is supposed turning clockwise around the origin. We notice 
that the three massive bodies are assumed to move in the same plane and that the model is not coherent in the sense that the 
primaries motions do not satisfy the Newton's equations. 

Let us take adimensional units as in the CR3BP and let mg = 328900.5614 be the mass of the Sun, ag be the distance between 
the Earth - Moon barycenter and the Sun, uj be the mean angular velocity of the Sun in synodical coordinates and 6>o be the 
value associated with the rotation of the Sun with respect to the Earth - Moon baycenter at t = 0. If 9 = ut, then the position of 
the Sun is described by 

xs = as cos (6» - 9 ), (5.1) 
Vs = -as sin (0 - ), 



and the equations of motion for the particle P can be written as 



(1 — u) . u ms „ „ ms 

x — 2y — x = — (x — (i) ^(x + 1 — fi) — (x — xs)—z cos (0 — 0q)—^-, 

rf r% rf, ai, 

y + 2x = y - V^fi-y - ±y - (y - ys )^f - sin (0 - 9 Q )^f, (5.2) 
r l r 2 r S a S 



i Z r 3' 
L r 1 



where fj, has the same meaning and value as the one introduced in Sec.[2]and r\ = [(a; — ^i) 2 + y 2 + z 2 ] 2 , r2 = [(x + 1 — /i) 2 + 
y 2 + z 2 \ 2 , rs = [(x — xs) 2 + {y — ys) 2 + z 2 ]? are the distances from P to Earth, Moon and Sun, respectively. 
We recall that this problem does not admit either first integrals or equilibrium points. 

We note that also in the case of a planet without a moon, it is still possible to apply the BR4BP by considering two planets and 
the Sun. For instance, we can assume the Sun and Mercury to move as in the CR3BP and Venus to move around their barycenter 
on a circular orbit lying on the same plane. For more details you can refer to [23 1 . 



5.1 Results 

Now, our objective is to clarify if the effect of the Sun can spoil relatively the heavier concentration of impact on the leading 
side of the Moon found previously. Hence, we apply the BR4BP equations of motion to the same initial conditions considered 
within the CR3BP framework. Also in this case, we are able to attribute to d^M some specific values, which account for the 
rate of recession of the Moon with respect to the Earth. We notice that as and u change accordingly to cLem^ as we assume the 
adimensional set of units defined in Sec. [3] The simulation is carried on as in Sec. [3j apart from the fact that now we have to 
explore the behaviour corresponding to different 8q (we take 5) and that we have to take care of impacts on the surface of the 
Earth. Finally, the maximum time span we allow to give birth to a lunar collision is of about 5 years. This choice is essentially 
due to the increasing computational effort. 
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Table 5.1 We show the percentages of Moon's and Earth's impact found, for different values of Earth - Moon distance d EM and initial 
phase of the Sun 9o. 



dEM.{km) 9o % Moon impacts % Earth impacts 



232400 


36° 


22.0 


2.7 


232400 


108° 


13.9 


3.3 


232400 


180° 


14.4 


3.0 


232400 


252° 


21.7 


2.6 



232400 324° 10.5 4.9 



270400 


36° 


20.0 


2.8 


270400 


108° 


10.1 


3.7 


270400 


180° 


10.5 


2.9 


270400 


252° 


20.1 


2.2 



270400 324° 6.8 4.2 



308400 


36° 


17.0 


2.4 


308400 


108° 


7.3 


3.7 


308400 


180° 


8.1 


2.1 


308400 


252° 


18.7 


2.0 



308400 324° 4.2 2.7 



384400 


36° 


13.3 


2.2 


384400 


108° 


3.2 


2.9 


384400 


180° 


3.9 


1.3 


384400 


252° 


14.8 


2.1 


384400 


324° 


1.2 


1.1 



Though we are aware that a deeper analysis using longer time intervales should be performed, significant results have already 
been obtained. They can be summarized as follows: 

(a) the percentage of impact depends on d^M ana< on the initial phase of the Sun, 8q: this is displayed in Tab. 15. 11 

(b) some trajectories collide with the Earth, the corresponding percentage is also shown in Tab. 15. II 
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Fig. 5.1 We display the density of impact (number of impacts per unit of area normalized with respect to the total number of impacts 
obtained) computed with the BR4BP equations of motion. The surface of the Moon is discretized in squares of 15° x 15°. The color bar 
indicates that the lighter the color of the square the greater the imp act d ensity. On the left, d = 36°; on the middle, B = 180° and on the 
right we display the distribution due to all the five values (see Tab. 15. 11 of Qo considered. On the top, d E M = 232400 km; on the bottom, 

d EM = 308400 km. 
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(c) looking to Tab. 15.11 it is also clear that there exist values of 6q more favorable to yield impacts with the Moon, while the 
Earth's percentage of impact has an almost constant trend; 

(d) it looks like the relative Earth - Moon and Earth - Moon - Sun distances, as well as the adimensional diameter of the Moon 
play a significant role in what concerns with the region of heavier lunar impact. In particular, the leading side collision 
concentration becomes more and more evident as cLem decreases. 

In Fig. 15. II we show the density of impact obtained when cLem = 232400 and cLem = 308400 km, respectively. For these 
plots, we consider 6q — 36°, 0q — 180° and the distribution deriving from all the values of 9q evaluated. 

6 A note on the computational effort 

As final remark, we note that all the simulation performed are quite expensive from a computational point of view, but that it 
was not an objective of this work to implement optimal programs. We use the UPC Applied Math cluster system, which consists 
in 26 Dell PowerEdge SC1425 servers, each with two Intel Xeon 3.2 GHz processors and 2 GB of RAM. In particular, we carry 
on parallel computations, using a node for each level of energy considered. In this way, the full computations for the CR3BP 
case take about 6 hours, for the BR4BP one about 1 day and the computations described in Sec.|4]about 1 hour and 50 minutes. 
In all the numerical integrations, we adopt a 7-8 Runge-Kutta-Fehlberg method. 

7 Conclusions and future developments 

The main purpose of this paper is to establish a relationship between low-energy trajectories in the Earth - Moon system and 
lunar impact craters. This is actually a quite wide and challenging topic, which involves knowledge related to mathematics, 
astronomy and geology. 

As primary goal, we define some tools which are effective in the determination of the dynamics pushing a massless particle 
under low-energy regimes. We exploit invariant objects within the Circular Restricted Three - Body Problem approximation, in 
particular transit trajectories lying inside the stable invariant manifold associated with the central invariant manifold of the Li 
equilibrium point. We implement a method that allows to reproduce the behaviour associated with the unstable component of 
any central orbit and does not need to distinguish between them. This fits with our investigation, because we are interested in 
minor bodies collisions that take advantage of the channels represented by the whole hyperbolic manifolds. To this purpose, we 
adopt well-known procedures to compute periodic Lyapunov orbits together with their hyperbolic invariant manifolds. 

With this approach, we perform extensive numerical simulations to determine both the lunar region of heavier impact and the 
sources of a potential uniform craters distribution. We also look for the influence of the Sun on these paths, by means of the Bicir- 
cular Restricted Four - Body Problem. Several outcomes can be highlighted, even if they have to be seen as patterns that require 
a more robust proof: further calculations with different dynamical and astronomical models are in progress. The investigation 
carried out is promising from many points of view, as it indicates future developments that are worth to be considered. 

Without the effect of the Sun, we get a confirmation that the neighbourhood of the apex of the surface of the Moon is the region 
where most collisions take place. We remark that the impact trajectories simulated reach the surface of the Moon with the lowest 
possible velocity: this point does not corrupt the apex concentration that other authors discovered without this restriction. The 
total time span considered (60 years) is sufficient to describe the general behaviour of the massless particles and no distinctions 
among different values for the Earth - Moon distance were observed. 

However, the gravitational force exerted by the Sun seems to blur the above phenomenon. Changing the ratio between the Earth 
- Moon - Sun distance and the Earth - Moon one, we notice different patterns. From our computations, it turns out that in more 
ancient epochs the low-energy lunar impacts were focused on the Moon leading side, but this is not true going further in time. 
Moreover, we get evidence that the position of the Sun with respect to the Earth - Moon barycenter affects the distribution of 
lunar impacts. These are the first aspects we plan to study with more detail. For instance, it would be deserving to integrate the 
BR4BP equations of motion for a longer interval of time and for a greater number of values of the initial phase of the Sun to 
understand the real nature of these numerical observations. 

On the other hand, we realize that small craters can also be generated by the impact of dust arising from more energetic 
collisions than the ones investigated here. Such phenomenon comes from the existence of periodic orbits that cross the surface 
of the Moon, that is, double collision orbits. In the future, we would like to see how they are transformed by the perturbation of 
the Sun. 

A natural step would be to add the gravitational attractions of other planets to see their consequences on the orbits simulated. 
This will be done by means of a Restricted n - Body Problem, using position and velocity of the primaries given by the JPL 
ephemerides (for instance the DE405 ones) and taking several initial epochs to compare the whole outcome. 

Moreover, we would like to link our methodology with real observational data, concerning either the existing lunar craters and 
the orbital parameters at a certain epoch of a given set of Near Earth Objects. This information would affect especially the way 
we generate the initial conditions corresponding to transit orbits. 
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Finally, to apply the same kind of analysis to the terrestrial planets would be of large interest. Starting from the CR3BP 
approximation, we mean to study the density of impact provided by Sun - planet low-enery orbits and then to add further 
gravitational effects, trying to figure out the orbital elements and also the regions in the phase space which more likely lead to 
collision. 
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Appendix: Knuth Shuffle Algorithm 

In this work, to obtain a sequence of random real numbers uniformly distributed between and 1 (U n ), we exploit the following linear congruential 
method. Let m = 2 31 - 1 = 2147483647, a = 7 s = 16807, q = 127773 and r = 2836, then 

T k = a(X k modq) - r[X k /q], 

( T k if T t >0 

aAi-mod m = < , 
\ T k + m if T k < 

= aX k mod m, 

C/fc+l = X k+1 /m. 

As shuffle algorithm, we mean a procedure which aims at reorder a given sequence (U n ) to improve its quality. We adopt this approach: 

1. we initialize an auxiliary sequence (Vb, V\, . . . , V p ) with the first p values of the X— sequence; 

2. we define Y = Xp+i; 

3. j = [pY/m\; 

4. Y = Vj- 

5. Vj = X p +i; 

6. the final output is represented by Y. 
For further details, please refer to 1161 . 
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